A New Template Family For The Detection Of Gravitational Waves From 
Comparable Mass Black Hole Binaries. 
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In order to improve the phasing of the comparable-mass waveform as we approach the last 
stable orbit for a system, various re-summation methods have been used to improve the standard 
post-Newtonian waveforms. In this work we present a new family of templates for the detection 
of gravitational waves from the inspiral of two comparable-mass black hole binaries. These new 
adiabatic templates are based on re-expressing the derivative of the binding energy and the 
gravitational wave flux functions in terms of shifted Chebyshev polynomials. The Chebyshev 
polynomials are a useful tool in numerical methods as they display the fastest convergence of any 
of the orthogonal polynomials. In this case they are also particularly useful as they eliminate one 
^ ' of the features that plagues the post-Newtonian expansion. The Chebyshev binding energy now 

has information at all post-Newtonian orders, compared to the post-Newtonian templates which 
only have information at full integer orders. In this work, we compare both the post-Newtonian 
and Chebyshev templates against a fiducially exact waveform. This waveform is constructed from 
a hybrid method of using the test-mass results combined with the mass dependent parts of the 
post-Newtonian expansions for the binding energy and flux functions. Our results show that the 
Chebyshev templates achieve extremely high fitting factors at all PN orders and provide excellent 
' parameter extraction. We also show that this new template family has a faster Cauchy convergence, 

gives a better prediction of the position of the Last Stable Orbit and in general recovers higher 
,1 ^ i Signal-to-Noise ratios than the post-Newtonian templates. 
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I. INTRODUCTION 



It is believed that stellar mass compact binaries consisting of double neutron stars (NS-NS), double black holes 
(BH-BH) or a mixed binary consisting of a neutron star and a black hole (BH-NS), are the primary targets for a 

■ direct first detection of gravitational waves (GW) by the ground-based interferometers, LIGO 1], VIRGO [2], GEO600 
O ■ 0, and TAMA Q]. It is also believed that the inspiral of comparable-mass Supermassive BHs during the merger of 
t**»« \ galaxies will be a major source of GWs for the proposed space-based detector LISA Under radiation reaction the 

■ orbit of a binary slowly decays, emitting a signal whose amplitude and frequency increases with time and is termed a 
"chirp" signal. While stellar models predict that there is a greater population of NS-NS binaries [E 0, IE IE EH , it is 
believed that BH-BH binaries will be the strongest candidates for detection since they can be seen within a greater 
volume, about two orders-of-magnitude greater than that for NS-NS binaries [EE]]- A quick calculation shows that 
the idealized Signal-to-Noise ratio (SNR) at 100 Mpc is ~ 2.5 for a NS-NS binary, ~ 5 for a BH-NS binary and ~ 12 
for a BH-BH binary, assuming the LIGO ground-based detector has a low frequency cutoff of 40 Hz. For LISA, it 
has been shown that we should, theoretically, be able to view the merger of Supermassive BHs out to cosmological 
distances of z ~ 10 [HE1. 

At present, the best proposed method to detect these sources is the method of matched filtering [l4j] . In this 
method, a set of waveforms or templates, that depend on a number of parameters of the source and its location and 
orientation relative to the detector, are created. These templates are then cross-correlated with the detector output 
weighted by the inverse of the noise spectral density. If a signal, whose parameters are close to one of the template 
waveforms, is present in the detector output then the cross-correlation builds up. In the case of a sufficiently strong 
signal, the correlation will be much larger than the root-mean-squared (RMS) correlation in the absence of any signal. 
The success of matched filtering depends on how well we understand the phase evolution of the waveform. A tiny 
instantaneous difference, as low as one part in 10 3 , between the phase of the true signal in the detector output and the 
search template, could lead to a cumulative difference of several radians as we integrate over several hundred to several 
thousand cycles. With the aim of improving the detection ratio for binary inspirals, there has been a considerable 
effort in accurately computing the dynamics of a compact binary and the emitted waveform. 

There has been a number of parallel efforts using different schemes to accurately describe the phase. Firstly, 
there is the post-Newtonian (PN) expansion of Einstein's equations which has been used to treat the dynamics 
of two bodies of comparable masses, with and without spin, in orbit around each other. This approximation is 
app licable when the velocities involved in the system are small but there is no restriction on the ratio of the masses 
HE EE E3, El El, HE El- Next, black hole perturbation theory has been used to compute the dynamics of a test 
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particle in orbit around a Schwarzschild or Kerr BH. Black hole perturbation theory does not make any assumptions 
on the velocit y o f the components, but is valid only in the limit when the mass of one of the bodies is much less than 
the other I22I I23L l2~il. l25l. \WL l27j . Finally, there has been some very promising progress made in the field of numerical 
relativity [281 ] - It is now possible to start numerical simulations a few orbits before merger, and evolve the black holes 
through the merger and into the ringdown phase. The goal is to take adiabatic inspiral templates and match them 
to the numerical waveforms as we approach the highly relativistic phase (29[. Hopefully, in time, this will provide us 
with a complete GW template, from the adiabatic regime of the inspiral through the merger and onto the ringdown. 

In this work, we are mostly interested in the PN approximation and the adiabatic inspiral phase. In particular we 
will focus on the inspiral of BH-BH binaries as these systems will coalesce in the most sensitive frequency band for 
the ground based detectors. The PN approximation is a perturbative method which expands the equations of motion, 
binding energy and GW flux as a power series in v/c, where v is the invariant velocity of the system and c is the 
speed of light. In the early adiabatic stage of an inspiral, the radiation reaction time-scale is much greater than the 
orbital time-scale. At this point in the binary evolution, the velocity of the bodies in the system is quite small. In 
fact for a BH-BH system of (10, 10) Mq, the velocity of the system as it enters the LIGO bandwidth is ~ 0.23 c. At 
this point, the two BHs are separated by a distance of R ~ 19 m, where m is the total mass of the system. This is 
quite close to the point at r < 10m (30j where we know the PN approximation begins to break down. 

It was shown (3ll . I32I HH, [H[ that templates based on re-summation methods, such as Pade approximation, have 
a faster convergence in modelling the gravitational waveform for test-mass systems. However, in this case an exact 
expression is known for the binding energy of the system, while the GW flux is known only in terms of a PN expansion 
in powers of v. When compared with numerical fluxes, the Pade approximations were clearly superior to the standard 
PN fluxes. The Pade based templates were then used to partially construct Effective One Body templates [35[ which 
went beyond the adiabatic approximation and modelled the waveform into the merg er phase. Other template families 
such as the more phenomenological BCV (Buonanno-Chen-Vallisneri) templates [36| and a new family called complete 
adiabatic [37j have also been proposed as possible successors to the PN template. 

In this paper, we build on the results of a previous work where we defined a new template family for the Schwarzschild 
test- mass case that was based on Shifted Chebyshev polynomials (SCPs). We were able to show that modelling the 
flux function with these orthogonal polynomials gave a more convergent and robust flux function which leads to lower 
errors in the estimation of the chirp mass. In this current work we explore the comparable-mass case. It is known that 
the Chebyshev polynomials are bound in the domain between ±1. In the previous work, we used the SCPs which we 
bound between zero and the velocity at the Last Stable Orbit (LSO). However, it is known in the literature that the 
smaller the domain we are approximating in, the faster the convergence of a Chebyshev series. We therefore derive 
the SCPs which are a function of the total mass of the system. This ensures that we minimize the domain of the 
SCPs for different comparable-mass systems to try and ensure the best results. Using the SCPs we propose a new 
representation of the binding energy and flux functions. 

A. Organization of the Paper. 

In Sec [H] we define the gravitational waveform in the restricted PN approximation for ground based detectors. We 
also define PN expansions for the derivative of the binding energy and the flux functions. In Sec lIIII we quickly review 
some aspects of Approximation theory and outline some of the details of the Chebyshev polynomials. We then go on to 
derive the analytic form for the SCPs in terms of the monomials in v, as well as analytic expressions for the monomials 
in terms of the SCPs. In Sec lIVI we define the new Chebyshev approximations to the binding energy and the GW flux 
function. We also define a fiducially exact binding energy and flux which is based on a hybrid scheme that uses aspects 
from both the test-mass and comparable-mass cases. Using these fiducial functions, we make a graphical comparison 
of the performance of both the PN and Chebyshev approximations against the respective fiducial functions. Sec [V] 
contains a more quantitative analysis of the new approximations. Using the newly defined functions we calculate the 
fitting factors of the Chebyshev and PN waveforms against a fiducial exact waveform. In Sec I VII we investigate the 
Cauchy convergence of this new template family. In Sees IVBl and I VIIH we look at how well the various approximations 
predict the position of the LSO, and how much signal-to-noise ratio they recover as a function of both total mass and 
distance when compared to the fiducial waveform. 

II. THE GRAVITATIONAL WAVEFORM. 

In the transverse-traceless gauge GWs are represented by the two polarizations h + and h x . The response h(t) 
of a detector to an incoming signal is given by the combination h(t) = h + F + + F x h x , where in the restricted 
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post-Newtonian approximation [381 ] the polarizations are given by 



, , , 4nm (1 + cos 2 i) ,,. 

M*) = 2 "W COS ^' W 

4 77 777. 

h x (t) = cosiu 2 (i) sin 0(f), (2) 

and the beam-pattern functions are given by 

F + (0,cf>,ip) = - (1 + cos 2 6) cos 20 cos 2?0 -cos sin 20 sin 2-0, (3) 

F x (6,(f>,ip) = - (1 + cos 2 0) cos 20 sin 2^ + cos sin 20 cos 2^. (4) 

By restricted PN approximation we retain only the dominant amplitude and ignore all higher harmonic corrections. 
In the above equations, m — mi + mi is the total mass of the system, r\ = 777,1777,2/m 2 is the reduced mass ratio, i is 
the inclination angle of the binary orbit, D is the distance to the source, v = (mfi) 1 / 3 = (irmfGw) 1 ^ 3 is the invariant 
PN velocity parameter related to the orbital frequency SI and GW frequency few, (8,4>) define the sky position of 
the source and ip is the polarization angle of the wave. As a single detector will be unable to disentangle the two 
individual polarizations, we can write the response in the form 

h(t) = -^-Av 2 (t) cos[0(f) + <fo], (5) 

where, for the short-duration signals we are considering in this study, the coefficient A = A(9, 0, i, ip) and phase factor 
0o = 0o (9, 0, i, ip) can be taken to be constant. 

The PN approximation and the quadrupole formula, when applied to an inspiralling binary system, give the rela- 
tivistic binding energy per unit mass E(v), and the GW flux F(v) as series expansions in the parameter v. Once we 
have these two functions we can use the energy balance argument 

m —jj— = -F(v) (6) 

to obtain the evolution of the phase of the GW. Integrating the energy balance equation and using 2irf = dep/dt, we 
obtain the set of parametrized integral equations 

t{v) = t +m dv—^-, (7) 
Jv F{v) 

0(«) = 0o + 2 / dvv 3 -f(, (8) 
Jv F{v) 

where E'(v) = dE(v)/dv, to and 0o is can be taken to be constants of integration at a particular reference time. 
Rather than solving the set of equations by numerical integration, it is much quicker to solve the following set of 
ordinary differential equations (ODEs) 

dv _ F(v) d(j> _ 2tj 3 

~db ~ ~mE'(v)' ~dt ~ ~m' ^ ' 

The above integral equations and ODEs for the phasing formula = 0(t) hold under the 'adiabatic inspiral' assump- 
tion. 

While for test-mass systems we have an exact expression for the binding energy and a PN expansion for the flux 
function, we have no such luxury in this case. For comparable mass bodies, the binding energy is represented by a 
PN series approximation to order 0(v 6 ), while the GW flux function is again represented by a series approximation 
to order 0(v 7 ). Thus we replace the functions (E'(v), F(v)) in Eqns (JH) with the PN power series approximations 
(E^(v),F Tn (v)). 

The PN energy derivative is defined by the power series expansion [39|, EfJ, \4M, [42, [43| 

6 

fc=0 



E'tAv) = - V vYe k v k (10) 
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and ei = 63 = e$ — 0. In the coefhcient ee we have used the fact that the constant term appearing in the literature 
has been determined to be A = —1987/3080 [4l|, [12, H, EH , which allows us to write the terms in square brackets in 
a more condensed form. 

The GW flux function, also defined in the form of a PN expansion, is given by [H, QJ, [HI [3 El S EU 
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where the Newtonian flux is given by 



F N (v) = ^r?v ia , 
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where 7 = 0.577... is Euler's constant. Once again, in coefficient Aq we have used the value for the constant 9 = 
— 11831/9240 [45| that appears in the literature to write the i] dependent term in a shortened form. We can see that 
the flux is not defined by a pure power series, but has a logarithmic term appearing at the 3-PN order. Terms like 
this can be responsible for poor convergence if not treated properly. 



III. THE SHIFTED CHEBYSHEV POLYNOMIALS AND THE VELOCITY MONOMIALS 

As we have seen, both the binding energy and the GW flux functions are given by power series representations to 
various orders. We can treat these expansions just like we would a standard Taylor series. We should caution the 
reader not to confuse the subscript T n with the symbol for the Chebyshev polynomial T n (x) that we will define in 
this section. While it is an unfortunate coincidence, we will continue to use both symbols as they are established in 
the literature. We have seen in other studies [U [3^, [H, Hi[ that in the case of a test-mass particle orbiting both 
Schwarzschild and Kerr black holes, the PN representation of the flux does not display the fastest convergence. Other 
re-summation methods such as Pade approximation offer a faster converging approximation to the GW flux function 
when compared with numerical results from Black Hole Perturbation Theory. It is these previous works that motivate 
us to try a different and more convergent approach to the modelling of the binding energy and flux functions. 
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The problem for expansions like a Taylor or a Maclaurin series is that they are based on the Weierstrass theorem, 
i.e. For any continuous function f(x) € C[a,b] and for any given e > 0, there exists a polynomial p n (x) for some 
sufficiently large n such that \f(x) — p n (x)\max < £■ So, as long as we can approximate a function with a series 
containing a sufficiently large number of terms, we can always minimize the error of our approximation. However, 
this is not always practical as the number of necessary terms may be too high (for example, it takes about 5000 terms 
in a Taylor expansion of arctan(x) to deliver a five significant figure accuracy at x equal to unity [46|). Moreover, 
we may be dealing with an approximation to some function where it is very difficult to calculate more than a few 
terms. For comparable mass systems this is a particular problem as we have only a three term expression for the 
binding energy and a seven term expansion for the GW flux function. We know from previous studies in the test- 
mass case that this number of terms may not be sufficient to accurately model the true binding energy and flux 
functions. A more promising possibility is based on the Chebyshev Alternation theorem for polynomials which states 
l For any continuous function f{x) € C[a,b], a unique minimax polynomial approximation p n (x) exists, and is uniquely 
characterized by the 'alternation property' (or 'equios dilation property') that there are (at least) n + 2 points in [a,b] 
at which \f{x) — p n {x)\ attains its maximum absolute value with alternating signs\ Thus, the reason the minimax 
polynomial is so sought after is that the Chebyshev Alternation theorem guarantees that there is only one such 
polynomial and it has the necessary condition of having an 'equal-ripple' error curve. Unfortunately, the minimax 
polynomial is usually extremely difficult, if not impossible, to find. A more promising and obtainable possibility is 
based on getting close to the minimax polynomial by using the family of Ultrasphcrical (or Gegenbauer) polynomials 
which are defined by 

P«(x) = C n (l - x*y a £; (1 - xT +a (-1 < a < oo) , (16) 

where C n is a constant. These polynomials are orthogonal over x G [—1, 1] with respect to the weight function 
(l — x 2 )° . A feature of the polynomials Pn \x) is that they have n distinct and real zeros and exhibit an oscillatory 
behavior in the interval [—1,1]. For the particular value of a = — 1/2 the amplitude of the oscillations remain constant 
throughout the interval and is conducive to finding an "equal-ripple" error curve, which is integral to the minimax 
polynomial. 

This value of a corresponds to the Chebyshev polynomials of the first kind T n (x) (hereafter Chebyshev polynomials). 
These polynomials are closely related to the minimax polynomial due to the fact that there are n + 1 points in [-1,1] 
where T n (x) attains a maximum value of unity with alternating signs [47j . It can be shown [48| that the Chebyshev 
polynomials exhibit the fastest convergence properties of all of the Ultraspherical polynomials. The n + 1 extrema 



are given by 

xk = cos [ — ) ft = 0, ..,n, (17) 

\n J 

with n zeros given by 

x k = cos P-* k=l,..,n. (18) 

The Chebyshev polynomials are formally defined by 

T n {x) = cos(n9) when x — cos(d). (19) 

From de Moivre's theorem, we know that cos (n9) is a polynomial of degree n in cos (6), e.g. 

cos (09) = 1, cos (6) = cos (9) , cos (29) = 2 cos 2 (9) - 1, , (20) 

which allows us to write the Chebyshev polynomials 

T (x) = 1, T x (x) = x, T 2 (x) = 2a; 2 - 1, (21) 



Therefore, a Chebyshev series in x corresponds to a Fourier series in 9. This close relation to a trigonometric function 
means these polynomials are extremely useful in approximating functions. The Chebyshev polynomials are orthogonal 

1/2 

polynomials with respect to the weight (l — x 2 ) according to 

(T i (x)\T j (x))= f dx Tii ^ Tji f = ( | i = Uo (22) 
J -l VI -x 2 [n i=j = 
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with initial conditions 



T (x) = l , T 1 {x)=x. 
We calculate the higher order Chebyshev polynomials using the recurrence relation 

T n (x) = 2xT n _ 1 {x)-T n _ 2 (x), 



(23) 
(24) 



given the above initial conditions. 

The Chebyshev polynomials are usually defined in the domain [-1,1]. For this particular problem, we require the 
Shifted Chebyshev Polynomials, T*(v), which are defined in the domain v <G [vini,vi so ] so that we have the maximum 
convergence possible. We can transform from the domain [-1,1] to another interval [a,b] using 



Now first defining the two variables 



2x -{a + b) 



X = Vini + Vu 



Vie [a,b],s £ [-1, 1]. 



£ = Vl so - Vi, 



(25) 



(26) 



where Vi n i = (irmfiow) 1 ^ 3 is the initial velocity of the body as it crosses the lower frequency cutoff of the detector and 
visa is the velocity at the last stable circular orbit, the shift has the form 



s = 



2v-x 



V V e [v in i,Vl so \. 



(27) 



We should point out here that both \ and £ are functions of the total mass of the binary through the initial velocity. 
This in turn infers that the Shifted Chebyshev polynomials (SCPs) are also functions of the total mass. Furthermore, 
while the PN approximations to the energy and flux functions are unchanged for all equal mass systems, regardless of 
total mass, this will not be the case for the Chebyshev approximations. For each equal mass system with a different 
total mass, the approximations to the energy and the flux will differ. We now write the SCPs in the form 



T*{v)=T n {s) = T n 



2v 



X 



and the recurrence relation as 



T*{v) 



2v-x 



l n-l 



(V) 



L n-2 



(V). 



(28) 



(29) 



While for computational purposes we will use the above recurrence relation to calculate the shifted polynomials, we 
also need to be able to express the polynomials in analytical form. Once we have the analytic expressions for T*{v) 
in terms of v, we can then invert the expressions to find an analytic format for the monomials of v in terms of T*(v). 
Therefore, using the fact that the SCPs are normalized according to 



T*(v) = 1 

T*(v) - rM2«-x] 

we can use the recurrence relation to find the other shifted polynomials to order 0(v 7 ), i.e. 
T*(v) = C 2 [(2x 2 - e) - &Xv - 8v 2 ] 

Ts(v) = C 3 [(3x£ 2 - 4 X 3 ) + 6 (4 X 2 - £ 2 ) v - A8 X v 2 + 32v 3 ] 

n(v) = r 4 [(£ 4 + 8x 4 - 8 X 2 e) + 32 ( X 2 e - 2 X 3 ) v + 32 (6 X 2 - £ 2 ) v 2 - 256 X v 3 + 218« 4 ] 



(30) 
(31) 

(32) 
(33) 
(34) 



Zb» = r 5 [(20 X 3 £ 2 - 16x 5 - 5 X £ 4 ) + (160 X 4 - 120 X 2 £ 2 + 10£ 4 ) v + (240 X £ 2 - 640 X 3 ) v 2 
+ (240x£ 2 - 640x 3 ) v 3 - 1280x« 4 + 512w 5 ] 



(35) 
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T;(v) = T 6 [(32x 6 - 48 X 4 £ 2 + 18 X 2 e 4 - f) + (384 X 3 f - 384 X 5 - 72 X £ 4 ) v + (l920 X 4 - 1152(x£) 2 + 72£ 4 ) v 2 

+ (l536x£ 2 - 5120x 3 ) v 3 + (7680x 2 - 768£ 2 ) v 4 ~ 6UA X v 5 + 2048v 6 ] (36) 

T;(v) = r 7 [(ll2x 5 C 2 ~ 64 X 7 - 56x 3 C 4 + 7x£ 6 ) + (896x 6 - 1120x 4 £ 2 + 336 X 2 £ 4 - 14£ 6 ) v 

+ (4480x 3 ? 2 - 5376x 5 - 672x£ 4 ) v 2 + (l7920x 4 - 8960(x£) 2 + 448£ 4 ) v 3 + (8960x£ 2 - 35840x 3 ) v 4 

+ (43008x 2 - 3584£ 2 ) v 5 - 28672x« 6 + 8192v 7 ] (37) 

Using the above expressions we can now write the monomials of v in terms of T* (v) as follows 

« = IbPo+Pi) (38) 

v 2 = (— + — \ T* + — T* + — T* (39) 



8 4 J u 2 1 8 

« 3 = ( 4xC 2 + 4) To + (lx 2 Z + I? + £rf a 2j + |U* (40) 



16^ 8 / u V 8 32" ) 1 16 /v " z 32 



b**fr)v*we***fa (42) 
" e - ( s + S- 4 ^ + ^ + da*') T » + (as** + S^ 3 + S*) T ' 

5S^* + 2li8 T » (43) 

- _ , X 7 35 6 _21_ 5 2 _L05_ 3 4 \ / 105 2 5 6 105 4 3 _35_. 7 \ r * 

128 + 2048 XC + 256 X 5 + 1024 X 5 J + ^ 1024 X § + 128 X C + 512 X * + 8192^ J 1 
21 5 2 _35_ 3 4 105 6 \ f 105_ 2 5 _35_ 43 21 7 \ 
256 X 5 + 256 X * + 496 X ^ J 2 + ^2048 X 5 + 512 X 5 + 8192^ J 3 



V 



V2048 1024^ ^ J 4 ^^2048^ * 8192" J 4096^" b 8192 

To express the PN expansions for both the energy derivative and the gravitational flux function in terms of shifted 
polynomials, we will substitute the above expressions for the monomials in v and collect all terms proportional to 
each polynomial. In Fig {]} we plot the SCPs for a (10, 10) M binary BH. We can see that the polynomials are 
bound between [-1,1] for v € [0.23, 0.408] which correspond to the initial and LSO velocities for this particular system 
assuming a lower frequency cutoff of 40 Hz and that the velocity at the LSO is given by the test-mass Schwarzschild 
value of viso = l/y/6. We can also see that SCPs have an equal amplitude oscillation exhibiting n + 1 alternating 
maximum and minimum values of ±1, as well as n zeros across the domain. 

IV. THE CHEBYSHEV APPROXIMATION GRAVITATIONAL WAVEFORM 

Our aim in this section is to take the set of ODEs given by Eqn ([§]) and write them in the form 

dv F Cn (v) dcj) 2v 3 



dt mE' c (vY dt 



(45) 
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FIG. 1: The 7th order Shifted Chebyshev Polynomials for a (10, 10) Mq BH-BH binary in the domain v G [0.23,0.408]. We 
can see that the polynomials have n + 1 alternating maximum and minimum values of ±1 and n zeros across the domain. 
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FIG. 2: A comparison of the convergence for the binding energy derivative coefficients of the PN approximation (left) and the 
shifted Chebyshev approximation (right) for the four BH-BH systems with masses (5, 5), (10, 10), (20, 10) and (20, 5) Mq. We 
can see that while the values of the PN coefficients oscillate wildly, the Chebyshev coefficients get smaller as we go to higher 
orders of approximation. We have plotted two systems with r\ = 0.25 to emphasize the fact that the Chebyshev coefficients 
£k ~ £k(m) are a function of the total mass of the system. 



where E' c (v) and Fc n (v) are more robust and convergent expansions of the binding energy and flux functions based 
onSCPs. 11 
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A. Modelling The Energy Function. 



We begin re- writing the monomials for v 2 ,v A and v 6 appearing in Eqn (fT0|) in terms of the SCPs. This allows us 
to define the Chebyshev approximation to the energy derivative function 



E' Cn ( V ) = - vv j2^n(v), 



where the new coefficients are defined by 



fc=0 



(46) 



So 
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ft = 6 6 ^ 5 
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Notice that the Chebyshev energy derivative is defined at all PN orders. The PN approximation suffered from the fact 
that for any odd PN order template, we had to use the previous order PN approximation for the energy derivative. 
This new expression allows us to use an equivalent order approximation to the flux function when we define waveforms 
at odd PN orders. We also remind the reader that the function E' c (v) = E' c (v; i], m) through the Shifted Chebyshev 
polynomials, whereas E' T (v) — E' T (v;rj) only. 

In Fig @ we plot the absolute value of the coefficients from both the PN and Chebyshev approximations for four 
BH-BH systems with masses (5, 5), (10, 10) and (20, 10) Mq to give values of r\ between 0.16 and 0.25. We can see from 
the left hand cell that the PN coefficients display no obvious convergence, with the values of the coefficients growing 
as we go to higher orders of approximation. On the other hand, the Chebyshev coefficients display a remarkable 
convergence with each successive coefficient smaller than the previous as we move to higher PN orders. This informs 
us that the Chebyshev series is convergent. We have shown in a previous work [49| that as the shifted polynomials have 
maximum absolute values of |T^(u)| < 1, the value of the coefficients £ n give an excellent estimate of the truncation 
error at each PN values for the Chebyshev series. 

In order to make some kind of comparison we have to define some fiducially exact energy derivative function. For 
this particular exercise we define a hybrid (Hoj energy function using the information we have from both the test-mass 
and comparable-mass cases. Using the 77 dependent terms from the PN energy we define 
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where the first term in the brackets is the exact expression for a test-particle in a Schwarzschild geometry. While we 
make the assumption that the comparable-mass case is a smooth 77 deformation of the test-mass case, this seems to 
be a plausible function to use as it returns the test-mass result in the case of 77 — > 0. 

In Fig ([3]) we plot the relative error between the PN and Chebyshev approximations and the fiducial energy function 
for a number of different BH-BH systems, i.e. 
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FIG. 3: Comparison of the relative errors for the PN (left) and Chebyshev (right) Energy Derivative functions against a 
fiducial exact energy function for four BH-BH systems with masses (5, 5), (10, 10), (20, 10) and (20, 5) Mq. Each dip in the 
approximation error curves correspond to a crossing of the exact energy function by the approximated flux. The initial velocity 
in each case is defined by a lower frequency cutoff of 40 Hz. 



We can see from the left hand column that in all cases the PN energy at 1-PN order has an expected high error. 
There is then a large jump to the 2-PN order and it is clear that the 3-PN energy gives a slightly better fit than the 
2-PN energy. From the images one could imagine that the PN approximation is converging to a particular solution. 
However, one only has to look at the PN approximation for the flux in both the test-mass and comparable-mass cases 
to realize that if it were possible to go to 4-PN order, it may be that the 4-PN approximation does worse than cither 
the 2 or 3-PN approximations. On the other hand, we observe a very remarkable convergence with the Chebyshev 
approximation. In most cases, the 1-PN order approximation is very close to the performance of the 3-PN Chebyshev 
approximation. For the low equal mass case, the approximation has quite a large oscillation at low velocity. However, 
this is increasingly minimized as we move to higher mass systems, both equal and unequal-equal. We can see that in 
most of the cases on the right hand side, it is virtually impossible to make out the various Chebyshev approximations 
from 1.5-PN order upwards. 



B. Modelling The Gravitational Wave Flux Function. 



In order to model the flux function as best as possible, we will borrow some steps from the work on Pade approx- 
imation [3l[. We can see from Eqn (fT2"|) that there is a null linear term in the PN expansion for the flux. We have 
confirmed that re- writing the expansion in terms of the SCPs automatically introduces a linear term in the expansion 
and that compared to the fiducial flux that we will define below, improves the approximation to the flux. However, 
we can do better. In the spirit of the Pade expansion, we first of all define the factored flux function 

/T n = (l-— W, ( 56 ) 

V v pole J 

where v po i e — l/-\/3 is the velocity at the light ring for a test-particle orbiting a Schwarzschild BH. In order to minimize 
the effect of the logarithmic term term appearing at 3-PN order, we factor out this term and normalize it according 
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PN Order 



FIG. 4: A comparison of the convergence for the flux coefficients of the PN approximation (left) and the shifted Chebyshev 
approximation (right) for the four BH-BH systems with masses (5, 5), (10, 10), (20, 10) and (20, 5)M@. We can see that the 
values of the PN coefficients again oscillate wildly and display an overall trend of growth. While not as smooth as in the case 
of the binding energy, the Chebyshev coefficients again get smaller as we go to higher orders of approximation. 
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where /o = cq and /& = A% — Ak-i/v po i e , for k — 1, . . . , 7. Now once again, swapping all of the monomials in v in 
the series expansion on the right hand side with their expansions in terms of the SCPs, we define the new Chebyshev 
flux function 
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where the new coefficients are defined by 
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FIG. 5: Comparison of the relative error for the PN (left) and Chebyshev (right) GW flux functions against a fiducial exact 
flux function for four BH-BH systems with masses (5, 5), (10, 10), (20, 10) and (20, 5) Mq. Each dip in the approximation error 
curve corresponds to a crossing of the exact energy function by the approximated flux. The initial velocity in each case is 
defined by a lower frequency cutoff of 40 Hz. 
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In a previous study [49j we showed that it was also prudent to expand the power after the log-term as a Chebyshev 
series. However, in that case we had a power series to order Oiv 11 ). In this study as the only information is given 
by the v e term, we found that there is no advantage to re-writing this term in terms of SCPs. Again we note that 
Fc n i v ) — Fc n (v;r),m) due to the dependence of the SCPs on v ini . In Fig (|4]) we again plot the absolute value of 
the flux coefficients for both the PN and Chebyshev approximations. We can see in this case that the PN coefficients 
(left) display large oscillations in their values as we go to higher PN orders. On the other hand, we once again see that 
the Chebyshev coefficients show convergent properties by getting smaller as we increase the order of approximation. 
We should point out that there is a slight oscillation in the values of the Chebyshev coefficients for the 77 = 0.16 case. 
The reason for this is not fully understood. What we have worked out is that this behaviour only exists for reduced 
mass ratios in the range 0.16 < i] < 0.19. 

As in the case of the energy function, we need to have some way of comparing both flux approximations. To this 
end, we use another hybrid scheme where we use the test-mass numerical flux from BH perturbation theory mixed 
with the r\ correction terms from PN theory. Thus the 3.5-PN fiducial exact flux is given by 
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where F^ n {v) is the Newtonian normalized numerical GW flux function (5lj |. 

In Fig iJSJ) we again plot the relative error between the PN and Chebyshev GW flux approximations and the fiducial 
flux function for the BH-BH systems of masses (5, 5), (10, 10), (20, 10) and (20, 5) M@. We can see from the cells in 
the left hand column a known trait of the PN expansion, that going to higher orders does not necessarily give a better 
approximation. In any event, it is clear that in all cases the 2.5-PN flux approximation is only marginally better than 
the 1-PN case. We should point out here that the 1.5 and 2-PN fluxes look much better than any other order as they 
approach the LSO. This is due to a coincidental crossing of the exact flux function as the LSO is reached. While the 3 
and 3.5-PN approximations are much better earlier on, it is clear that their relative errors are greater than the 1.5 and 
2-PN orders at the LSO. However, we should point out that for the (20, 5) Mq case, the 3-PN approximation has a 
crossing as the LSO approaches, thus giving a very small error in the most important region. The PN approximations 
can give an extremely good fit in a very narrow range as the approximated flux crosses the exact flux. However, it is 
clear from the plot that where this happens for each PN order depends on the total mass of the system. 

The Chebyshev fluxes, on the other hand, display a remarkable convergence. In all cases we see oscillating error 
curves indicating that the Chebyshev approximations cross the exact flux many times. While the 1-PN case has the 
largest amplitude of error oscillation, all other orders of approximation are incredibly consistent. We can see that for 
the equal mass cases, the Chebyshev fluxes at all orders have a smaller error than PN fluxes at the 3 and 3.5-PN 
order. It is clear to see from the slopes of the error curves, that the Chebyshev approximation attempts to find an 
equal amplitude error curve by leaving the error float in some parts in order to improve the performance elsewhere. 
While it does not fully succeed, it is pretty much impossible to make out any difference between the Chebyshev 
approximations at PN orders of greater than one. While the Chebyshev approximations do not have the small range 
accuracy of the PN approximations, overall the Chebyshev approximations display a higher stability. 



V. FITTING FACTORS. 

We define the overlap between two waveforms h(t) and s(t) as the inner product of the normalized waveforms 
denoted by 

O = , {kls) , (68) 

where the scalar product is defined by 

{h[s) = 2 [ ' Mf) [ h w r w +h *w s w]- (69) 

Here, the * denotes complex conjugate and h(f), s(f) are the Fourier transforms of h(t), s(t). The above scalar 
product is weighted by the inverse one-sided noise power spectral density (PSD) of the detector Sh(f)- For initial 
LIGO, the design study PSD Q is given by [32j 

Sh(f) = 9 x 10~ 46 [0.52 + 0. 16x" 4 - 52 + 0.32a; 2 ] Hz" 1 , (70) 

where x = f / fk, and //. = 150 Hz is the "knee- frequency" of the detector. We take the PSD to be infinite below 
the lower frequency cutoff of /i ow = 40 Hz. For testing the performance of a particular template, a more appropriate 
quantity compared to the scalar product is the fitting factor FF. Defining each template as a function of the 
intrinsic parameters A'' = {to, rj, to, (j)o}, the fitting factor is defined as the maximum overlap obtained by varying the 
parameters of the template relative to the signal we are trying to detect. 

FF = maxO (A p ) . (71) 

While we expect the comparable mass waveforms to have LSO frequencies different to the test-mass case, we have 
verified that we can take the fi so from the test-mass case as the upper limit to the scalar product integral without 
any great change of results. For this study we generate our signal by inserting the equations for the fiducially exact 
binding energy and flux functions given by Equations (|54[) and (|67|) into the set of ODEs that describe the velocity 
and phase evolution of the waveform, i.e. hx(v) — h{E x (v), Fx(v)). 
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FIG. 6: The percentage errors in total mass m, reduced mass ratio r/ and fitting factors for three equal mass systems of 
(5, 5), (10, 10) and (20, 20) Mq for both PN and Chebyshev waveforms when compared to the fiducially exact hybrid waveform. 




FIG. 7: The percentage errors in total mass m, reduced mass ratio r\ and fitting factors for three unequal mass systems of 
(10, 5), (20, 10) and (20, 5) Mq for both PN and Chebyshev waveforms when compared to the fiducially exact hybrid waveform. 



In Figs ((6]) and (fTj> we plot the percentage errors in the recovered total mass m and reduced mass ratio rj, as well 
as the recovered fitting factors for six different comparable mass systems. For this study we have chosen three equal 
mass systems of (5, 5), (10, 10) and (20, 20) Mq, and three unequal mass systems of (10, 5), (20, 10) and (20, 5) Mq. 
Firstly, focusing on the PN waveforms. We can see that in all cases the 1-PN template performs very badly against 
the test waveform, never achieving very high fitting factors and always returning extremely high errors in the total 
mass (almost 90% in the case of the (20.5) Mq system.). From 1.5-PN order onwards, we can always achieve fitting 
factors of greater than 0.95. However the figures display the known oscillatory nature of the PN expansion. We can 
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see that the 2-PN is better than the 1.5-PN, but then the 2.5-PN order template is worse than the 2-PN template 
etc. In all cases however, the 3 and 3.5-PN order templates return acceptable answers. We should bring attention 
to the fact that for the (5, 5) Af© system, the 2-PN template outperforms the Chebyshev template in the extraction 
of the parameters, with a comparable fitting factor to the Chebyshev templates. However, this behaviour does not 
repeat itself in any of the other cases we examined. 

On the other hand, the Chebyshev templates display a remarkable convergence with all templates always achieving 
a fitting factor of greater than 0.99 . In terms of parameter estimation, it is clear, in all cases, that the error at 1-PN 
order is only slightly greater than the Chebyshev errors at 3.5-PN order. This is a consequence of the Chebyshev 
series attempting to minimize the maximum error and find a minimax solution. While there is some oscillation in the 
error curves, it is quite small in comparison to the PN templates. We can see from the figures, that compared to the 
PN templates, the Chebyshev templates always achieve comparable or smaller errors in the estimation of m and r\. 



VI. CAUCHY CONVERGENCE. 



If a sequence {x n } converges, the terms get closer and closer to the limit of the sequence as the order of the 
approximation increases. However, the Cauchy criterion demands that rather than a limit, the terms get closer to 
each other. To this extent, we define our Cauchy convergence as 



C = (h n \h n+ i) 



(72) 



In calculating the Cauchy convergence, the masses of both templates are kept the same as we only need to maximize 
over the extrinsic parameters t and 4>q . Maximization over t is achieved by simply computing the correlation of the 
template with the data in the frequency domain followed by the inverse Fourier transform. This yields the correlation 
of the signal with the data for all time-lags. To maximize over (f> we use the definition of the "Best Possible Overlap" 
when individually maximizing over the phases of two separate templates [3lj . This is given by 
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The Cauchy convergence has been calculated many times before for PN templates (See for example (31(), so we 
will not repeat the exercise here. Suffice to say that the PN templates display the usual oscillatory behaviour which 
shows that going from one PN order to the next does not necessarily result in a better template. In Table (|TJ we list 
the Cauchy convergence for the Chebyshev templates for the six test cases we analysed in the previous sections. The 
parameter n denotes the level of approximation, e.g. n = 2 corresponds to 1-PN or v 2 . We can see from the tables 
that the Chebyshev templates are incredibly convergent with all templates achieving overlaps of greater than 0.9 with 
the successive template. In fact it is only in the (5, 5) Mq case that the template fails to achieve overlaps of > 0.97 
with each other. 



VII. POSITION OF THE LSO. 



In the previous sections of this work we took the position of the LSO to be at R = 6 m thus giving the dimensionless 
velocity of vi so — Here we investigate how close both the PN and Chebyshev approximations come to correctly 

predicting the position of the LSO, defined as E'(v) = 0. For both the 'exact' and PN expansions this quantity is 
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n 


(5, 5)M 


(10, 1O)M 


(20, 2O)M 


(10,5)M Q 


(20, 10)M o 


(20, 5)M Q 


2 


0.9497 


0.9841 


0.99991 


0.9703 


0.9998 


0.99987 


3 


0.99969 


0.9999 


1.0 


0.99953 


0.999996 


0.99978 


4 


0.9998 


1.0 


1.0 


0.99949 


1.0 


0.99998 


5 


0.9999 


1.0 


1.0 


0.99999 


1.0 


1.0 


6 


1.0 


1.0 


1.0 


1.0 


1.0 


1.0 



TABLE I: The Cauchy convergence {h n \h n +i) of the Chebyshev templates for three equal mass, and three unequal mass systems 
assuming a LIGO noise curve and a lower frequency cutoff of 40 Hz. 

E'(v;rj) = 0, while for the Chebyshev approximation it is E'(v;T),m) = giving different values for different total 
masses. 

In Fig ([5]) we plot the positions of the predicted LSOs against the fiducial value for four of the BH-BH systems that 
we used previously. In each cell, the circle denotes the value of the exact LSO position in units of R/m. We have 
placed the exact position of the LSO at 3-PN to signify the order of r\ corrections involved. The squares denote the 
PN approximation values and the triangles denote the Chebyshev approximations. 

In the top two cells we use the same values for both the fiducial and PN predicted LSOs. These cells correspond 
to the rj = 0.25 case which does not change for the PN approximation, but does for the Chebyshev approximation. 
We can see that in all four cases, the PN approximation at 1-PN gives a very bad prediction of the position of the 
LSO. This value is improved somewhat at 2-PN order, with the best value coming in all cases from the 3-PN order 
PN approximation. 

Just like other re-summation methods, there is no new information gained from restructuring the series approxima- 
tion. Thus, there are no obvious tricks we can use to ensure that the location of the LSO is improved at the highest 
PN order. Therefore, it is no surprise that the PN and Chebyshev templates predict the same location at 3-PN order. 
However, we can see from the image that the prediction for the LSOs position is greatly improved at lower PN orders 
with the Chebyshev templates. The flat profile of the error distribution is consistent with the effort of the Chebyshev 
approximation to minimize the maximum error. 

VIII. SIGNAL-TO-NOISE RATIOS. 

The final aspect we will investigate is a comparison of recovered signal-to-noise ratios (SNRs) using the fiducially 
exact template as our signal. Using the definition of the scalar product, we can write the signal-to-noise ratio (SNR) 
obtained by searching for a signal hx (t) with a template h n (t) as 

p = (K \hx) , 7g j 
\J (h n \ h n ) 

For this investigation we place the exact waveform and template at a distance of 100 Mpc. We also assume that the 
parameters of the template match the parameters of the signal exactly. In Fig we plot the recovered SNRs at 
various PN orders for equal mass systems. The results scale as y/^rj for unequal mass binaries. We can see from the 
plots that at 1, 1.5 and 2.5-PN orders the Chebyshev templates recover a much greater SNR than the PN templates 
at virtually all masses. At the 2-PN order the Chebyshev templates outperform the PN templates at low masses and 
are then comparable at the higher masses. While for the 3 and 3.5-PN orders, the PN templates outperform the 
Chebyshev templates at very low masses, but are then comparable at higher masses. 

IX. CONCLUSION. 

In this work we have presented a new family of templates to detect GWs from the inspiral of comparable mass 
BH-BH binaries. The templates are based on shifted Chebyshev polynomials which are bound between the initial and 
final velocities of the system in question. In general, the Chebyshev polynomials display the fastest convergence of 
all the orthogonal polynomials. Using the shifted polynomials, we have defined a new re-summed binding energy and 
GW flux function which we demonstrated to be graphically more convergent than the PN functions when compared 
with a fiducially exact energy and flux. Using a fiducially exact template constructed from the exact energy and 
flux, we were able to show that the Chebyshev templates achieve very high fitting factors and excellent parameter 
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FIG. 8: Comparison of the position of the LSO for the PN and Chebyshev energy functions against the minimum of the fiducial 
exact energy function for the four BH-BH systems with masses (5, 5), (10, 10), (20, 10) and (20, 5) Mq. The circles denote the 
LSO for the exact fiducial energy, the squares represent the PN approximation and the triangles represent the Chebyshev 
approximation. We can see that the Chebyshev approximation is more convergent and gives a better prediction of the location 
of the LSO at virtually all PN orders. 



extraction as compared to the PN templates. We were also able to show that in an effort to get close to a minimax 
solution, the 1-PN Chebyshev template is almost as good as the 3.5-PN PN template. Further displaying the strength 
of these new templates, we were also able to show that they provided more accurate measurements for the position 
of the LSO, have an incredibly fast Cauchy convergence, and in most cases achieve a higher recovered SNR. 

We should point out that the Chebyshev templates outperformed the PN templates when compared to a fiducially 
exact waveform constructed from a hybrid combination of the test-mass binding energy and a numerical flux function 
from BH perturbation theory. These functions were then combined with the mass dependent parts of the PN expansion 
for the energy and flux functions for comparable-mass bodies. While is seems sensible, as we recover the test-mass 
result in the limit rj — ► 0, it is still only an approximation of a true possible waveform. However, we are confident in the 
abilities of this new template family. We plan to follow up this investigation with a comparison against comparable- 
mass waveforms from numerical relativity. This should provide us with a more concrete critique of the ability of these 
new Chebyshev templates to detect GWs and parametrize their sources. 

For now, we believe that this new family of comparable-mass templates will be a valuable addition to the search 
for BH-BH binaries with both ground and space-based detectors. 
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